Introducción a R como GIS
Introducción a R como GIS
- Resumen
- Asignar el directorio de trabajo
- Trabajando con vectores desde R
- Descargar datos espaciales desde R
- Trabajando con rasters desde R
- La belleza de lo interactivo
- Un ejemplo de programación funcional con
purrr
## Loading required package: sp
## Checking rgeos availability: TRUE
Resumen
En este tutorial vermos como utilizar R como un Sistema de Información Geográfica (SIG). R utiliza un conjunto de librerÃas para procesar objetos espaciales; las liberÃas que vamos a utilizar en este tutorial son:
- sp: Classes and methods for spatial data; the classes document where the spatial location information resides, for 2D or 3D data. Utility functions are provided, e.g. for plotting data as maps, spatial selection, as well as methods for retrieving coordinates, for subsetting, print, summary, etc.
- rgdal: Provides bindings to Frank Warmerdam’s Geospatial Data Abstraction Library (GDAL) (>= 1.6.3) and access to projection/transformation operations from the PROJ.4 library. The GDAL and PROJ.4 libraries are external to the package, and, when installing the package from source, must be correctly installed first. Both GDAL raster and OGR vector map data can be imported into R, and GDAL raster data and OGR vector data exported. Use is made of classes defined in the sp package. Windows and Mac Intel OS X binaries (including GDAL, PROJ.4 and Expat) are provided on CRAN.
- raster: Reading, writing, manipulating, analyzing and modeling of gridded spatial data. The package implements basic and high-level functions. Processing of very large files is supported.
- rgeos: Interface to Geometry Engine - Open Source (GEOS) using the C API for topology operations on geometries. The GEOS library is external to the package, and, when installing the package from source, must be correctly installed first. Windows and Mac Intel OS X binaries are provided on CRAN.
- leaflet: Create and customize interactive maps using the ‘Leaflet’ JavaScript library and the ‘htmlwidgets’ package. These maps can be used directly from the R console, from ‘RStudio’, in Shiny apps and R Markdown documents.
- rasterVis: Methods for enhanced visualization and interaction with raster data. It implements visualization methods for quantitative data and categorical data, both for univariate and multivariate rasters. It also provides methods to display spatiotemporal rasters, and vector fields. See the website for examples.
En general trabajaremos sobre tareas especÃficas y vermos como ejecutarlas desde R. Finalmente se asume que los estudiantes ya tiene conocimientos básicos de SIG como:
Vector
- Qué es un vector
- Cuantos tipos de vectores hay?
- Puntos
- LÃneas
- PolÃgonos
- Conocer la estructura interna de un vector
- Número de vectores
- Tabla de atributos
Raster
- qué es un raster?
- Qué es un extent?
- Cuáles son sus dimensiones?
- Saber qué es una proyección?
- Saber si un raster está proyectado o no?
Asignar el directorio de trabajo
Trabajando con vectores desde R
Como se mencionó existen tre tipos de vectores espaciales; estos son puntos, lÃneas y polÃgonos. R utiliza un conjunto de librerÃas para procesar objetos espaciales. Las librerÃas que utilizamos para leer y hacer operaciones sobre vectores son: rgdal y sp
Cómo leo un vector desde R.
El comando para leer un shapefile es readOGR; este se encuentra en la librerÃa rgdal. Los argumentos del comando son:
dsn: Es el folder donde tenemos nuestro archivo vector (shapefile)layer: El nombre de la capa que queremos leer
## rgdal: version: 1.3-6, (SVN revision 773)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
## Linking to sp version: 1.3-1
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/Luis/Dropbox/CursoNichos2019/destdv1gw", layer: "destdv1gw"
## with 382 features
## It has 8 fields
## Integer64 fields read as strings: COV_ COV_ID
Grafiquemos el mapa
Explorando los atributos del vector
Uno de los comandos de más útilies para explorar los elementos que conforman a un objerto, es el comando str.
o bien para darnos solo una idea de su estrutura podemos escribir el nombre del objeto
Los objetos espaciales en R son de la clase S4; para extraer los atributos de estos objetos (indexar) se utiliza el sÃmbolo de **@**. Veamos como extraer la tabla de atributos del mapa map_vec
| AREA | PERIMETER | COV_ | COV_ID | ENTIDAD | CAPITAL | RASGO_GEOG | CVE_EDO | |
|---|---|---|---|---|---|---|---|---|
| 0 | 6.7179224 | 18.9211577 | 2 | 1 | BAJA CALIFORNIA | Mexicali | NA | 02 |
| 1 | 16.6952450 | 30.6009734 | 3 | 5 | SONORA | Hermosillo | NA | 26 |
| 2 | 0.0000470 | 0.0318774 | 4 | 2 | BAJA CALIFORNIA | NA | ISLA | 02 |
| 3 | 0.0000177 | 0.0175657 | 5 | 3 | BAJA CALIFORNIA | NA | ISLA | 02 |
| 4 | 0.0001608 | 0.0679551 | 6 | 4 | BAJA CALIFORNIA | NA | ISLA | 02 |
| 5 | 0.0127180 | 0.5160367 | 7 | 8 | BAJA CALIFORNIA | NA | ISLA | 02 |
Extraer un polÃgono del shapefile
puebla <- map_vec[map_vec@data$ENTIDAD=="PUEBLA",]
puebla <- map_vec[map_vec@data$ENTIDAD %in% "PUEBLA",]
plot(puebla)Guardar un shapefile
Ahora guardamos como un shapefile separado, el polÃgono
Filtrar los puntos que caen dentro de un polÃgono
Generemos unos puntos aleatorios definidos en el extent del mapa; el comando nombre_objeto@bbox nos muestra el extent de nuestro shapefile.
## min max
## x -118.36603 -86.71074
## y 14.53401 32.71877
# Coordenadas aleatorias
rand_x <- runif(n =300, map_vec@bbox[1,1],max = map_vec@bbox[1,2])
rand_y <- runif(n =300, map_vec@bbox[2,1],max = map_vec@bbox[2,2])
coords_rand <- data.frame(longitude=rand_x,latitude=rand_y)Transformamar puntos a coordenadas espaciales
Puntos en el poligono
data_poly_all <- sp::over( sp_df_rand , map_vec , fn = NULL)
en_poligono_index <- which(!is.na(data_poly_all[,1]))
p_en_poligono <- sp_df_rand[en_poligono_index ,]
plot(map_vec)
plot(p_en_poligono,add=T,col="blue")Extraer los atributos en los puntos del poligono
data_en_poli_extract <- data_poly_all[en_poligono_index,]
coords_en_poli <- sp_df_rand@coords[en_poligono_index,]
data_en_poli_extract_df <- cbind(coords_en_poli,data_en_poli_extract)
sp_DF_en_poli_extract <- SpatialPointsDataFrame(data_en_poli_extract_df[,c("longitude","latitude")],
data_en_poli_extract_df)
#writeOGR(sp_DF_en_poli_extract,
# dsn = map_dir,layer = "puntos_aleatorios_mexico",
# driver = "ESRI Shapefile")| longitude | latitude | AREA | PERIMETER | COV_ | COV_ID | ENTIDAD | CAPITAL | RASGO_GEOG | CVE_EDO | |
|---|---|---|---|---|---|---|---|---|---|---|
| 2 | -110.09291 | 29.04229 | 16.695245 | 30.60097 | 3 | 5 | SONORA | Hermosillo | NA | 26 |
| 10 | -97.82715 | 19.84400 | 2.939128 | 18.83949 | 298 | 313 | PUEBLA | Herorica Puebla de Zaragoza | NA | 21 |
| 14 | -97.09354 | 19.36795 | 6.062359 | 33.15909 | 255 | 269 | VERACRUZ DE IGNACIO DE LA LLAVE | Xalapa-Enriquez | NA | 30 |
| 20 | -103.97467 | 27.88340 | 22.882945 | 29.38276 | 10 | 10 | CHIHUAHUA | Chihuahua | NA | 08 |
| 22 | -98.87711 | 23.35904 | 6.917870 | 30.18376 | 71 | 75 | TAMAULIPAS | Ciudad Victoria | NA | 28 |
| 24 | -93.17818 | 16.89131 | 6.214071 | 17.55246 | 371 | 389 | CHIAPAS | Tuxtla Gutierrez | NA | 07 |
Reproyectar un polÃgono
Con R podemos hacer reproyección de coordenadas, podemos pasar de coordenadas geográficas (logitude, latitud) a coordenadas utm (metrosal nivel del mar).
Veamo como proyectar nuesto mapa de Puebla a coordenadas utm
puebla_reproj <- spTransform(puebla,CRS("+proj=utm +zone=10+datum=WGS84"))
plot(puebla_reproj,axes=TRUE)Regresemos a coordenadas geográficas
puebla_reproj_geo <- spTransform(puebla_reproj,CRS("+proj=longlat +ellps=WGS84"))
plot(puebla_reproj_geo,axes=TRUE)Descargar datos espaciales desde R
Ahora veremos como descargar datos espaciales desde R; para ello utlizaremos la función getData del paquete raster. La función permite descargar datos de diferentes bases de datos como:
LÃmites territoriales
- GADM: Datos de lÃmites territoriales de paÃses
Datos climáticos
- worldclim: Datos climáticos
Datos Digitales de Elevación
Trabajando con rasters desde R
La librerÃa estrella para trabjar con capas raster es raster. Esta es un wrapper que se comunica con funciones de alto nivel escritas en el lenguaje C para procesar objetos raster. Es importante notar que los raster no son nada más que arreglos bidimensionales de tipo numérico o categórico.
Al igual que los archivos de tipo vector tienen un extent, sin embargo tienen otros atributos como dimensión y resolución espacial. Entre mayor sea la resolución del raster, mayor será su dimensión (mayor número de filas y columnas).
Cómo leer un raster
El comando para leer 1 solo raster es raster
Cómo leer un conjunto (stack) de capas
A veces (casi siempre) será necesario leer un conjunto de capas que tienen la misma resolución y extent (i.e. WorldClim layers). El comando stack de raster permite crear stack de capas; para ello es necesario indicar con en un vector los path o rutas a la carpeta que contiene nuestras capas.
Veamos como leer o almacenar en la memoria de R el conjunto de capas de WorldClim layers
Primero en guardamos en una variable los path a nustras capas
Usamos el comando stack
Graficamos las primer 5 capas
notemos que utilice la misma sintaxis que utilizamos para indexa listas…
Hacer un crop
La función crop utiliza el extent de un shapfile para hacer recortar un raster o un stack. Veamos un recorte de las capas de WorldClim con el poligono de puebla.
Hacer una mascara
Al igual que crop, mask hace un recorte pero conserva la forma del poligono.
Guardar un raster
El comando para gurardar un raster es wirteRaster y está en el paquete raster
Guardar un stack completo
Mostraremos como guardar un stack completo en una carpeta
Convertir a puntos un stack
| x | y | bio1 | bio10 | bio11 | bio12 | bio13 | bio14 | bio15 | bio16 | bio17 | bio18 | bio19 | bio2 | bio3 | bio4 | bio5 | bio6 | bio7 | bio8 | bio9 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| -97.75000 | 20.75000 | 240 | 274 | 192 | 1265 | 202 | 37 | 56 | 561 | 125 | 522 | 133 | 106 | 53 | 3329 | 332 | 134 | 198 | 270 | 202 |
| -97.75000 | 20.58333 | 237 | 270 | 189 | 1737 | 282 | 49 | 64 | 811 | 150 | 784 | 165 | 109 | 54 | 3287 | 330 | 130 | 200 | 266 | 199 |
| -97.91667 | 20.41667 | 222 | 254 | 178 | 2278 | 421 | 48 | 77 | 1179 | 151 | 905 | 161 | 119 | 56 | 3054 | 321 | 111 | 210 | 248 | 188 |
| -98.25000 | 20.25000 | 148 | 171 | 122 | 1219 | 246 | 23 | 77 | 630 | 81 | 298 | 115 | 134 | 65 | 1923 | 249 | 45 | 204 | 157 | 133 |
| -98.08333 | 20.25000 | 176 | 203 | 141 | 2142 | 402 | 42 | 79 | 1144 | 144 | 480 | 149 | 128 | 60 | 2506 | 277 | 67 | 210 | 194 | 151 |
| -97.91667 | 20.25000 | 199 | 228 | 159 | 2587 | 480 | 55 | 77 | 1359 | 180 | 985 | 190 | 121 | 58 | 2794 | 298 | 91 | 207 | 221 | 169 |
Extraer valores de un raster o un stack
Podemos extraer los valores de un stack utilizando coordenadas.
coord_rand_extract <- data.frame(extract(bios_wc,coords_rand))
coord_rand_extract <- na.omit(coord_rand_extract)
knitr::kable(head(coord_rand_extract))| bio1 | bio10 | bio11 | bio12 | bio13 | bio14 | bio15 | bio16 | bio17 | bio18 | bio19 | bio2 | bio3 | bio4 | bio5 | bio6 | bio7 | bio8 | bio9 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2 | 231 | 293 | 161 | 465 | 144 | 3 | 110 | 313 | 19 | 279 | 80 | 158 | 51 | 5217 | 381 | 75 | 306 | 288 | 225 |
| 9 | 167 | 256 | 74 | 271 | 54 | 6 | 71 | 135 | 20 | 122 | 26 | 160 | 45 | 7076 | 338 | -13 | 351 | 245 | 125 |
| 10 | 171 | 196 | 137 | 1832 | 372 | 42 | 70 | 884 | 144 | 413 | 158 | 124 | 61 | 2358 | 267 | 66 | 201 | 187 | 147 |
| 13 | 221 | 290 | 139 | 639 | 119 | 23 | 50 | 242 | 92 | 182 | 98 | 130 | 43 | 5880 | 362 | 64 | 298 | 265 | 150 |
| 14 | 149 | 169 | 124 | 1721 | 321 | 41 | 76 | 872 | 127 | 460 | 160 | 120 | 66 | 1800 | 238 | 57 | 181 | 160 | 133 |
| 15 | 175 | 259 | 82 | 683 | 88 | 28 | 37 | 235 | 96 | 212 | 96 | 145 | 42 | 6877 | 333 | -5 | 338 | 225 | 82 |
Tabla de correlaciones…
| bio1 | bio10 | bio11 | bio12 | bio13 | bio14 | bio15 | bio16 | bio17 | bio18 | bio19 | bio2 | bio3 | bio4 | bio5 | bio6 | bio7 | bio8 | bio9 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| bio1 | 1.00 | 0.65 | 0.87 | 0.35 | 0.41 | 0.07 | 0.08 | 0.39 | 0.07 | 0.33 | 0.14 | -0.52 | 0.28 | -0.44 | 0.45 | 0.88 | -0.60 | 0.56 | 0.68 |
| bio10 | 0.65 | 1.00 | 0.20 | -0.02 | -0.14 | 0.21 | -0.33 | -0.15 | 0.21 | -0.10 | 0.19 | -0.34 | -0.48 | 0.40 | 0.88 | 0.26 | 0.16 | 0.49 | 0.34 |
| bio11 | 0.87 | 0.20 | 1.00 | 0.44 | 0.60 | -0.06 | 0.33 | 0.60 | -0.07 | 0.47 | 0.04 | -0.43 | 0.68 | -0.82 | 0.03 | 0.98 | -0.87 | 0.40 | 0.69 |
| bio12 | 0.35 | -0.02 | 0.44 | 1.00 | 0.90 | 0.71 | -0.32 | 0.91 | 0.70 | 0.82 | 0.77 | -0.66 | 0.22 | -0.43 | -0.30 | 0.55 | -0.62 | -0.09 | 0.40 |
| bio13 | 0.41 | -0.14 | 0.60 | 0.90 | 1.00 | 0.37 | 0.08 | 0.99 | 0.36 | 0.83 | 0.45 | -0.50 | 0.49 | -0.65 | -0.35 | 0.64 | -0.73 | 0.04 | 0.43 |
| bio14 | 0.07 | 0.21 | -0.06 | 0.71 | 0.37 | 1.00 | -0.81 | 0.38 | 0.99 | 0.51 | 0.94 | -0.64 | -0.39 | 0.17 | -0.08 | 0.10 | -0.12 | -0.27 | 0.18 |
| bio15 | 0.08 | -0.33 | 0.33 | -0.32 | 0.08 | -0.81 | 1.00 | 0.08 | -0.82 | -0.06 | -0.70 | 0.51 | 0.67 | -0.50 | -0.11 | 0.17 | -0.20 | 0.25 | 0.07 |
| bio16 | 0.39 | -0.15 | 0.60 | 0.91 | 0.99 | 0.38 | 0.08 | 1.00 | 0.37 | 0.82 | 0.47 | -0.48 | 0.49 | -0.65 | -0.35 | 0.63 | -0.72 | 0.02 | 0.43 |
| bio17 | 0.07 | 0.21 | -0.07 | 0.70 | 0.36 | 0.99 | -0.82 | 0.37 | 1.00 | 0.51 | 0.95 | -0.64 | -0.41 | 0.18 | -0.08 | 0.09 | -0.11 | -0.25 | 0.18 |
| bio18 | 0.33 | -0.10 | 0.47 | 0.82 | 0.83 | 0.51 | -0.06 | 0.82 | 0.51 | 1.00 | 0.59 | -0.53 | 0.32 | -0.49 | -0.37 | 0.52 | -0.63 | 0.03 | 0.36 |
| bio19 | 0.14 | 0.19 | 0.04 | 0.77 | 0.45 | 0.94 | -0.70 | 0.47 | 0.95 | 0.59 | 1.00 | -0.62 | -0.28 | 0.07 | -0.09 | 0.19 | -0.21 | -0.27 | 0.33 |
| bio2 | -0.52 | -0.34 | -0.43 | -0.66 | -0.50 | -0.64 | 0.51 | -0.48 | -0.64 | -0.53 | -0.62 | 1.00 | 0.13 | 0.21 | 0.12 | -0.61 | 0.60 | -0.16 | -0.36 |
| bio3 | 0.28 | -0.48 | 0.68 | 0.22 | 0.49 | -0.39 | 0.67 | 0.49 | -0.41 | 0.32 | -0.28 | 0.13 | 1.00 | -0.92 | -0.42 | 0.57 | -0.70 | 0.03 | 0.35 |
| bio4 | -0.44 | 0.40 | -0.82 | -0.43 | -0.65 | 0.17 | -0.50 | -0.65 | 0.18 | -0.49 | 0.07 | 0.21 | -0.92 | 1.00 | 0.49 | -0.77 | 0.90 | -0.08 | -0.45 |
| bio5 | 0.45 | 0.88 | 0.03 | -0.30 | -0.35 | -0.08 | -0.11 | -0.35 | -0.08 | -0.37 | -0.09 | 0.12 | -0.42 | 0.49 | 1.00 | 0.00 | 0.44 | 0.43 | 0.20 |
| bio6 | 0.88 | 0.26 | 0.98 | 0.55 | 0.64 | 0.10 | 0.17 | 0.63 | 0.09 | 0.52 | 0.19 | -0.61 | 0.57 | -0.77 | 0.00 | 1.00 | -0.90 | 0.38 | 0.71 |
| bio7 | -0.60 | 0.16 | -0.87 | -0.62 | -0.73 | -0.12 | -0.20 | -0.72 | -0.11 | -0.63 | -0.21 | 0.60 | -0.70 | 0.90 | 0.44 | -0.90 | 1.00 | -0.15 | -0.55 |
| bio8 | 0.56 | 0.49 | 0.40 | -0.09 | 0.04 | -0.27 | 0.25 | 0.02 | -0.25 | 0.03 | -0.27 | -0.16 | 0.03 | -0.08 | 0.43 | 0.38 | -0.15 | 1.00 | 0.09 |
| bio9 | 0.68 | 0.34 | 0.69 | 0.40 | 0.43 | 0.18 | 0.07 | 0.43 | 0.18 | 0.36 | 0.33 | -0.36 | 0.35 | -0.45 | 0.20 | 0.71 | -0.55 | 0.09 | 1.00 |
Reproyectar un raster
Al igual que con los vectores, podemos reproyectar a los rasters. Lo anterior se puede hacer mediante el comando projectRaster de raster
puebla_bio1 <- puebla_ras[[1]]
puebla_bio1_reproj <- projectRaster(puebla_bio1,
crs="+proj=utm +zone=10+datum=WGS84")
plot(puebla_bio1_reproj)
plot(puebla_reproj,add=T)Regresemos a geográficas
puebla_bio1_geo <- projectRaster(puebla_bio1,
crs="+proj=longlat +ellps=WGS84")
plot(puebla_bio1_geo)
plot(puebla,add=T)La belleza de lo interactivo
En este ejemplo veremos como limpiar datos de precias de una especie con un mapa interactivo
Leemos los datos
| name | longitude | latitude | prov | date | key |
|---|---|---|---|---|---|
| Puma concolor | -80.68902 | 25.57658 | gbif | 2019-01-02 | 1986571407 |
| Puma concolor | -108.56825 | 33.12471 | gbif | 2019-01-04 | 1986581111 |
| Puma concolor | -116.57024 | 32.74562 | gbif | 2019-01-05 | 1986595849 |
| Puma concolor | -93.47792 | 16.05619 | gbif | 2018-01-08 | 1831091833 |
| Puma concolor | -120.73842 | 35.32380 | gbif | 2018-01-30 | 1837077255 |
| Puma concolor | -120.25755 | 39.43723 | gbif | 2018-01-24 | 1807289765 |
- Quitamos datos duplicados
library(ntbox)
puma_occsDF_L1 <- ntbox::clean_dup(puma_occsDF,
longitude = "longitude",
latitude = "latitude",
threshold = 0) Visualizamos los datos
puma_occsDF_L1$rcID<- as.character(1:nrow(puma_occsDF_L1))
map %>% leaflet::addMarkers(lng = puma_occsDF_L1$longitude,
lat = puma_occsDF_L1$latitude,
popup = puma_occsDF_L1$rcID)Limpiamos datos
d_raros <- c(791,1367,966)
puma_occsDF_L2 <- puma_occsDF_L1[-d_raros ,]
map %>% leaflet::addMarkers(lng = puma_occsDF_L2$longitude,
lat = puma_occsDF_L2$latitude,
popup = puma_occsDF_L2$rcID)Escribimos nuestra tabla de datos
Descarga de datos de presencia spocc
library(spocc)
jaguar_occs <- occ("Panthera onca",
from = "gbif",
limit = 1500,
has_coords = T)
jaguar_occsDF <- occ2df(jaguar_occs)
map %>% leaflet::addCircleMarkers(lng = jaguar_occsDF$longitude,
lat = jaguar_occsDF$latitude)Un poco más profesional
# install.packages("wicket")
library(wicket)
mex <- wrld_simpl[wrld_simpl$NAME %in% c("Mexico"),]
raster::extent(mex)## class : Extent
## xmin : -118.4042
## xmax : -86.7014
## ymin : 14.55055
## ymax : 32.71846
jaguar_mexExtent <- occ(query = 'Panthera onca',
geometry =c(-118.40417,14.55055,
-86.70140,32.71846))
jaguar_mexExtentDF <- occ2df(jaguar_mexExtent)
map %>% leaflet::addCircleMarkers(lng = jaguar_mexExtentDF$longitude,
lat = jaguar_mexExtentDF$latitude)Aún más pro
Un ejemplo de programación funcional con purrr
- Leer la base de datos de reptiles y anfibios
- Partir el
data.framepor especie. - Crear polÃgnos mÃnimos convexos de los registros.
- Cortar las capas de
WorldClimcon el polÃgono. - Escribir los rasters (capas de modelación).
- Buscar los datos de gbif usando polÃgonos.
- Limpiar los duplicados espaciales.
- Regresar el raster y los puntos de gbif.
- Partir el
##
## Attaching package: 'magrittr'
## The following object is masked from 'package:raster':
##
## extract
Definimos la función para crear y escribir polÃgonos convexos en el ambiente global
library(rgdal)
library(raster)
library(spocc)
# Objeto (funcion) en el ambiente global
create_convex <- function(sp_db,sp_name,long,lat,dir_sv){
ch <- chull(sp_db[,c(long,lat)])
coords <- sp_db[c(ch, ch[1]),
c(long,lat)]
sp_poly <- SpatialPolygons(list(Polygons(list(Polygon(coords)),
ID=1)))
sp_poly_df <- SpatialPolygonsDataFrame(sp_poly,
data=data.frame(ID=1))
writeOGR(sp_poly_df, dir_sv,
layer=sp_name,
driver="ESRI Shapefile",
overwrite_layer = TRUE,delete_dsn=TRUE)
return(sp_poly_df)
}Crando la función
mi_funcion <- function(sp_db,sp_name,long,lat,dir_sv,env_layers,n_occs=500){
dir_sv <- as.character(dir_sv)
sp_name <- as.character(sp_name)
if(!dir.exists(dir_sv)) dir.create(dir_sv)
mi_ch <- create_convex(sp_db,sp_name,long,lat,dir_sv)
sp_Layerscortadas <- crop(env_layers, mi_ch)
sp_Layerscortadas <- mask(sp_Layerscortadas, mi_ch)
lapply(names(sp_Layerscortadas), function(x){
writeRaster(sp_Layerscortadas[[x]], paste0(dir_sv,"/",
x,".asc"),overwrite=TRUE)})
sp_occs <- occ(sp_name,
from = "gbif",
limit = n_occs,
has_coords = T,
geometry = wicket::sp_convert(mi_ch))
sp_occsDF <- occ2df(sp_occs)
sp_occsLimpio <- ntbox::clean_dup(sp_occsDF,
longitude = "longitude",
latitude = "latitude",
threshold = raster::res(sp_Layerscortadas)[1])
listaRes <- list(sp_occs=sp_occsLimpio,
layers= sp_Layerscortadas,
poligono=mi_ch)
return(listaRes)
}Veamos con una especie
d_sps_DF <- d_spsL[[1]]
sp_ejemplo <- mi_funcion(sp_db = d_sps_DF,
sp_name = d_sps_DF$name[1],
long = "longitude",
lat = "latitude",
dir_sv = d_sps_DF$name[1],
env_layers = wc10,n_occs = 500)Grafiquemos
plot(sp_ejemplo$layers[[1]])
plot(sp_ejemplo$poligono,add=T)
points(sp_ejemplo$sp_occs[,c("longitude","latitude")],
col="blue",pch=20)Ahora hagámoslo con algunas (ustedes pueden con todas las sps)
spsLista <- d_spsL[1:3]
spsResL <- spsLista %>% purrr::map(~mi_funcion(sp_db = .x,
sp_name = .x$name[1],
long = "longitude",
lat = "latitude",
dir_sv = .x$name[1],
env_layers = wc10,
n_occs = 500))Exploremos la lista
plot(spsResL$`Anaxyrus cognatus`$layers[[1]])
plot(spsResL$`Anaxyrus cognatus`$poligono,add=T)
points(spsResL$`Anaxyrus cognatus`$sp_occs[,c("longitude","latitude")],
col="red",pch=20)